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Abstract 

We consider a stocliastic Susceptible- Exposed-Infected- Recovered (SEIR) epidemiological model. 
Through the use of a normal form coordinate transform, we are able to analytically derive the 
stochastic center manifold along with the associated, reduced set of stochastic evolution equations. 
The transformation correctly projects both the dynamics and the noise onto the center manifold. 
Therefore, the solution of this reduced stochastic dynamical system yields excellent agreement, 
both in amplitude and phase, with the solution of the original stochastic system for a temporal 
scale that is orders of magnitude longer than the typical relaxation time. This new method allows 
for improved time series prediction of the number of infectious cases when modeling the spread of 
disease in a population. Numerical solutions of the fluctuations of the SEIR model are considered 
in the infinite population limit using a Langevin equation approach, as well as in a finite population 
simulated as a Markov process. 
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The reduction of high-dimensional, stochastic systems is an important and 
fundamental problem in nonlinear dynamical systems. In this article, we present 
a general theory of stochastic model reduction which is based on a normal form 
coordinate transform. This nonlinear, stochastic projection allows for the deter- 
ministic and stochastic dynamics to interact correctly on the lower-dimensional 
manifold so that the dynamics predicted by the reduced, stochastic system agree 
well with the dynamics predicted by the original, high-dimensional, stochastic 
system. Although the method may be applied to any physical or biological 
system with well-separated time scales, here we apply the method to an epi- 
demiological model. We show that when compared to the original, stochastic 
epidemic model, the reduced model properly captures the initial and recurrent 
disease outbreaks, both in amplitude and phase. This long-term accuracy of the 
reduced model allows for the application of effective disease control where phase 
differences between outbreak times and vaccine controls are important. Addi- 
tionally, in practice, one can only measure the number of infectious individuals 
in a population. Our method allows one to predict the number of unobserved 
exposed individuals based on the observed number of infectious individuals. 



I. INTRODUCTION 



The interaction between deterministic and stochastic effects in population dynamics has 
played, and continues to play, an important role in the modeling of infectious diseases. The 
mechanistic modeling side of population dynamics is well-known and establishedii^. These 
models typically are assumed to be useful for infinitely large, homogeneous populations, and 
arise from the mean field analysis of probabilistic models. On the other hand, when one con- 
siders finite populations, random interactions give rise to internal noise effects, which may 
introduce new dynamics. Stochastic effects are quite prominent in finite populations, which 
can range from ecological dynamics^ to childhood epidemics in cities^i^. For homogeneous 
populations with seasonal forcing, noise also comes into play in the prediction of large out- 
breaks^i^. Specifically, external random perturbations change the probabilistic prediction 
of epidemic outbreaks as well as its control^. 

When geometric structure is applied to the population, the interactions are modeled as a 
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networki^iii. Many types of static networks which support epidemics have been considered. 
Some examples include small world networks^^, hierarchical networks^^, and transportation 
networks of patch models^. In addition, the fluctuation of epidemics on adaptive networks, 
where the wiring between nodes changes in response to the node information, has been 
examined^. In adaptive network models, even the mean field can be high-dimensional, 
since nodes and links evolve in time and must be approximated as an additional set of 
ordinary differential equations. 

Another aspect of epidemic models which is often of interest involves the inclusion of a 
time delay. The delay term makes the analysis significantly more complicated. However, 
it is possible to approximate the delay by creating a cascade consisting of a large number 
of compartments^. For example, one could simulate the delay associated with a disease 
exposure time with several hundred "exposed" compartments. 

These model examples are just a few of the types of very high-dimensional models that 
are currently of interest. As a result of the high- dimensionality, there is much computation 
involved, and the analysis is quite difficult. In particular, real-time computation is not 
currently possible. However, there are usually many time scales that are well-separated (due 
typically to a large range in order of magnitude of the parameters) when considering such 
high-dimensional problems. In the presence of well-separated time scales, a model reduction 
method is needed to examine the dynamics on a lower- dimensional space. It is known that 
deterministic model reduction methods may not work well in the stochastic realm, which 
includes epidemic models^. The purpose of this article is to examine a method of nonlinear, 
stochastic projection so that the deterministic and stochastic dynamics interact correctly on 
the lower-dimensional manifold and predict correctly the dynamics when compared to the 
full system. Because the noise affects the timing of outbreaks, it is essential to produce a 
low-dimensional system which captures the correct timing of the outbreaks as well as the 
amplitude and phase of any recurrent behavior. 

We will demonstrate that our stochastic model reduction method properly captures the 
initial disease outbreak and continues to accurately predict the outbreaks for time scales 
which are orders of magnitude longer than the typical relaxation time. Furthermore, in 
practice, real disease data includes only the number of infectious individuals. Our method 
allows us to predict the number of unobserved exposed individuals based on the observed 
number of infectious individuals. 
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For stochastic model reduction, there exist several potential methods for general prob- 
lems. For a system with certain spectral requirements, the existence of a stochastic center 



manifold was proven in Ref. 
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Non-rigorous stochastic normal 
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form analysis (which leads 
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Rigorous theoretical 



to the stochastic center manifold) was performed in Refs. 
analysis of normal form coordinate transformations for stochastic center manifold reduction 
was developed in Refs. |23|]24 Later, another method of stochastic normal form reduction 
was developed^^, in which any anticipatory convolutions (integrals into the future of the noise 
processes) that appeared in the slow modes were removed. Since this latter analysis makes 
the construction of the stochastic normal form coordinate transform more transparent, we 
use this method to derive the reduced stochastic center manifold equation. 

Figure [T] shows a schematic demonstrating our approach to the problem. We consider a 
high-dimensional system along with its corresponding reduced low-dimensional system. In 
this article, two types of low-dimensional system are discussed: a reduced system based on 
deterministic center manifold analysis and a reduced system based on a stochastic normal 
form coordinate transform. Regardless of the type of low- dimensional system being con- 
sidered, a common noise is injected into both the high- dimensional and low-dimensional 
models, and an analysis of the solutions found using the high and low-dimensional systems 
is performed. 




High-Dimensional 
System 



Associated 
Low— Dimensional 
System 



FIG. 1: Schematic demonstrating the injection of a common noise into both the high-dimensional 
system and its associated low-dimensional system. 



In this article, as a first study of a high-dimensional system, we consider the Susceptible- 
Exposed-Infected-Recovered (SEIR) epidemiological model with stochastic forcing. As pre- 
viously mentioned, we could easily consider a SEIR-type model where the exposed class was 
modeled using hundreds of compartments. Since the analysis is similar, we consider the 
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simpler standard SEIR model to demonstrate the power of our method. Section [TTl provides 
a complete description of this model. Section IIIII describes how to transform the determin- 
istic SEIR system to a new system that satisfies the spectral requirements needed to apply 
the center manifold theory. After the theory is used to find the evolution equations that 
describe the dynamics on the center manifold, we show in Sec. IIVI how the reduced model 
that is found using the deterministic result incorrectly projects the noise onto the center 
manifold. Section |V] demonstrates the use of a stochastic normal form coordinate transform 
to correctly project the noise onto the stochastic center manifold. A discussion section and 
the conclusions are contained respectively in Sec. |VT] and Sec. IVIII 



II. THE SEIR MODEL FOR EPIDEMICS 



We begin by describing the stochastic version of the SEIR model found in Ref. |26|. We 
assume that a given population may be divided into the following four classes which evolve 
in time: 

1. Susceptible class, s{t), consists of those individuals who may contract the disease. 

2. Exposed class, e(t), consists of those individuals who have been infected by the disease 
but are not yet infectious. 

3. Infectious class, i(t), consists of those individuals who are capable of transmitting the 
disease to susceptible individuals. 

4. Recovered class, r{t), consists of those individuals who are immune to the disease. 

Furthermore, we assume that the total population size, denoted as A^, is constant and can 
be normalized to S{t) + E{t) + I{t) + R{t) = 1, where S{t) = s{t)/N, E{t) = e{t)/N, 
I{t) = i{t)/N, and R{t) = r{t)/N. Therefore, the population class variables S, E, I, and R 
represent fractions of the total population. The governing equations for the stochastic SEIR 
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model are 



I{t) = aE{t) - (7 + /i)J(t) + a303(t), 



E{t) = pI{t)S{t) - (a + fi)E{t) + (T^Mt), 



S{t) = fi- pI{t)S{t) - fiS{t) + (7101 (t), 



(la) 
(lb) 
(Ic) 



R{t) = jl{t) - fiR{t) + a^Mt), 



(Id) 



where cXi is the standard deviation of the noise intensity Di = crf/2. Each of the noise terms, 
0i, describes a stochastic, Gaussian white force that is characterized by the correlation 
functions 



Additionally, fi represents a constant birth and death rate, (3 is the contact rate, a is the 
rate of infection, so that 1/a is the mean latency period, and 7 is the rate of recovery, so 
that 1/7 is the mean infectious period. Although the contact rate jS could be given by a 
time-dependent function (e.g. due to seasonal fluctuations), for simplicity, we assume P to be 
constant. Throughout this article, we use the following parameter values: fi = 0.02(year)~^, 
/3 = 1575.0(year)~^, a = l/0.0279(year)~^, and 7 = l/O.O^year)"-"^. Disease parameters 
correspond to typical measles values^^i^^. Note that any other biologically meaningful pa- 
rameters may be used as long as the basic reproductive rate Ro = aP/[{a + /i)(7 + /i)] > 1. 
The interpretation of Rq is the number of secondary cases produced by a single infectious 
individual in a population of susceptibles in one infectious period. 

As a first approximation of stochastic effects, we have considered additive noise. This 
type of noise may result from migration into and away from the population being consid- 
ered^S. Since it is difficult to estimate fluctuating migration rates^E, it is appropriate to treat 
migration as an arbitrary external noise source. Also, fluctuations in the birth rate manifest 
itself as additive noise. Furthermore, as we are not interested in extinction events in this 
article, it is not necessary to use multiplicative noise. In general, for the problem considered 
here, it is possible that a rare event in the tail of the noise distribution may cause one or 
more of the S, E, and I components of the solution to become negative. In this article, we 
will always assume that the noise is sufficiently small so that a solution remains positive for 
a long enough time to gather sufficient statistics. Even though it is difficult to accurately 




(2a) 




(2b) 
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estimate the appropriate noise level from real data, our choices of noise intensity lie within 



29|. The case for multiplicative noise will be 



the huge confidence intervals computed in Ref. 
considered in a separate paper. 

Although S + E + I + R= lin the deterministic system, one should note that the 
dynamics of the stochastic SEIR system will not necessarily have all of the components sum 
to unity. However, since the noise has zero mean, the total population will remain close 
to unity on average. Therefore, we assume that the dynamics are sufficiently described by 
Eqs. ( p^ - lITcl) . It should be noted that even if E(t) + I(t) = for some t, the noise allows 
for the reemergence of the epidemic. 



III. DETERMINISTIC CENTER MANIFOLD ANALYSIS 

One way to reduce the dimension of a system of equations is through the use of determin- 
istic center manifold theory. In general, a nonlinear vector field can be transformed so that 
the linear part (Jacobian) of the vector field has block diagonal form where the first matrix 
block has eigenvalues with positive real part, the second matrix block has eigenvalues with 
negative real part, and the third matrix block has eigenvalues with zero real part^Si^i. These 
blocks are associated respectively with the unstable eigenspace, the stable eigenspace, and 
the center eigenspace. If we suppose that there are no eigenvalues with positive real part, 
then orbits will rapidly decay to the center eigenspace. 

In order to make use of the center manifold theory, we must transform Eqs. (fTa|) - (fTc|) to a 
new system of equations that has the necessary spectral structure. The theory will allow us 
to find an invariant center manifold passing through the fixed point to which we can restrict 
the transformed system. Details regarding the transformation can be found in Sec. IIII A[ 
and the computation of the center manifold can be found in Sec. IIIIBI 



A. Transformation of the SEIR model 



Our analysis begins by considering the governing equations for the stochastic SEIR model 
given by Eqs. (ITal) - (fTc|) . We neglect the crj0j(t) terms in Eqs. (ITa|) - (ITc|) so that we are 
considering the deterministic SEIR system. This deterministic system has two fixed points. 
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The first fixed point is 

(^e,i?e,/e) = (1,0,0), (3) 

and corresponds to a disease free, or extinct, equilibrium state. The second fixed point 
corresponds to an endemic state and is given by 

/ (7 + /i) (a + /i) /i /i (7 + /i) fia n\ 

\ pa a + fi ap (7 + /i) (a + /i) p J 

To ease the analysis, we define a new set of variables, S, E, and J, as S{t) = S{t) — So, 
E{t) = E{t) — Eo, and I(t) = I{t) — Jq. These new variables are substituted into Eqs. fllap - 

Then, treating /i as a small parameter, we rescale time by letting t = fir. We may then 
introduce the following rescaled parameters: a = ao/fi and 7 = 7o//i, where ao and 70 
are 0{1). The inclusion of the parameter /i as a new state variable means that the terms 
in our rescaled system which contain /i are now nonlinear terms. Furthermore, the system 
is augmented with the auxiliary equation ^ = 0. The addition of this auxiliary equation 
contributes an extra simple zero eigenvalue to the system and adds one new center direction 
that has trivial dynamics. The shifted and rescaled, augmented system of equations is 

dS (ao + /i^)(7o + ^ aofJ'^P ^ x 

— = -PfiIb ■ — r-- ■ — (5aj 

dr ao {ao + n^){'jo + fi^) 

f =MS + + + 7 + " + ^'^^V ^'^^ S-{ao + ^^^)E, (5b) 

dr ao ("0 + /^ )(7o + /^ ) 

dl 



^ =aoE - (70 + (5c) 



where the endemic fixed point is now located at the origin. 

The Jacobian of Eqs. (I5al)-( l5dl) is computed to zeroth-order in /i and is evaluated at 
the origin. Ignoring the fi components, the Jacobian has only two linearly independent 
eigenvectors. Therefore, the Jacobian is not diagonalizable. However, it is possible to 
transform Eqs. fl5al) - fl5cl) to a block diagonal form with the eigenvalue structure that is 
needed to use center manifold theory. We use a transformation matrix, P, consisting of the 
two linearly independent eigenvectors of the Jacobian along with a third vector chosen to be 
linearly independent. There are many choices for this third vector; our choice is predicated 
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on keeping the vector as simple as possible. This transformation matrix is given as 



00+70 
70 







"0+70 g ]^ 
70 



(6) 



Using the fact that {S, E, I)^ = P ■ (f/, V, W)'^, then the transformation matrix leads to the 
following definition of new variables, U, V, and W: 



ao + 7o 



= ^ + 



7o 



"0 + 7o 
W = I + E. 



-E. 



(7a) 
(7b) 
(7c) 



The application of the transformation matrix to Eqs. ( I5al) -(l5cl) leads to the transformed 
evolution equations given by 



^ = -aoU+ (^0^ - "0^) _ (70 + /^') («o + /i') [(«o + lo) U + IqW] 
dr ao + 7o ao (ao + 7o) 



"0 + 7o 



ao7o 



(7o + A^^) ("0 + /^^) 



iU + V) 



(8a) 



dV ^^^^ _ fi^ (70 1/ - apU) _ (70 + /i^) (qq + /i^) [(qq + 70) + -fpW] 
dr ao + 7o 70 (oq + 70) 



7o ("0 + 7o) 



7oiy + (tto + 7o) U 



At «o7o 



(70 + At^) ("0 + A^^ 



(8b) 



(ir 



«ot/- (70 + Ai') (f/ + l^) + 



(70 + /x2) (ao + A^') [(ao + 7o) ^ + 7oW^] 



ao7o 



A^ "070 



;,2y + AZIZ + {ao + ^o)U 

7o V (70 + 1^'') (ao + A^^) 



{U + V) 



(8c) 



=0. 



(8d) 
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B. Center manifold equation 



The Jacobian of Eqs. (!8a|) - (18dp to zeroth-order in /i and evaluated at the origin is 



-("0 + 7o) 



(ao+7o) 



«070 



(ao+7o) 



(9) 



which shows that Eqs. fl8aj) - fl8d|) may be rewritten in the form 

— = Ax + f(x,y,/i), 



dy 

dr 



By + g(x,y,/i), 



(10) 

(11) 
(12) 



where x = (U), y = {V, W), A is a constant matrix with eigenvalues that have negative real 
parts, B is a constant matrix with eigenvalues that have zero real parts, and f and g are 
nonlinear functions in x, y and /i. In particular. 



A 



-("0 + 7o) 



B 











CtQ-fO 







(13) 



Therefore, the system will rapidly collapse onto a lower-dimensional manifold given by 
center manifold theory^^. Furthermore, we know that the center manifold is given by 



U = h{V,W,fi), 



(14) 



where h is an unknown function. 
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Substitution of Eq. (fT^ into Eq. (l8a|) leads to the following center manifold condition: 
dh{V,W,fi)dV , dh{V,W,fi)dW , fi^joV-aoh{V,W,fi)] 



dV dr dW dr u v , . ao + 7o 

(70 + /i') (gp + /i') [(qq + 70) h{V, W, /i) + 70 TV] 
"0 ("0 + 7o) 

^oW^ + («o + 70) h{V, W, /i) + J^^f^-—] {h{V, W, fi) + V). (15) 



«0 + 7o V ' ' (70 + yU^) ("0 + yU^) , 

In general, it is not possible to solve the center manifold condition for the unknown function, 
h(y, W, jj). Therefore, we perform the following Taylor series expansion of /i(V, VF, /x) in V , 
W, and /i: 

h{V, W, /i) =ho + h^V + hjW + h^^l + h22V^ + /issV^Vr + /l33Vr' + 

/i^2/iV^ + /i^S/"^ + V/"' + • • • ' (16) 

where /iq, /i2, ^3, ^/x, • • • are unknown coefficients that are found by substituting the Taylor 
series expansion into the center manifold condition and equating terms of the same order. 
By carrying out this procedure using a second-order Taylor series expansion of /i, the center 
manifold equation is 

where e = \{y, W, /i)| so that e provides a count of the number of V, W, and /i factors in any 
one term. Substitution of Eq. (fT7|) into Eqs. (l8bl) and (IHcI) leads to the following reduced 
system of evolution equations which describe the djTiamics on the center manifold: 

dV /i^7o^aoVr fi'^aoW 'JofJ^'^V 



dr {ao + 7o)^ (ao + lof "0 + 7o 

{,o + ,^)aoW _ Pa^, ,Hao + ,o) \ lo^W 

«o + 70 (ao + 7o) V (7o + /i^) (^o + ^J^^) J \ (qq + 70)' 

dW fi^W , fi^W 2t 



.■2+— + 

dr {ao + 7o) ao + 7o 

^W+ ^ /^M«o + 7o) \(v-^^^]. (18b) 



"0 + 70 V (70 + At^) (ao + At^) / V (ao + 7o) 
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IV. INCORRECT PROJECTION OF THE NOISE ONTO THE STOCHASTIC 
CENTER MANIFOLD 



A. Transformation of the stochastic SEIR model 

We now return to the stochastic SEIR system given by Eqs. (fTa|) -( !Tc|l . The shift of the 
fixed point to the origin will not have any effect on the noise terms, so that the stochastic 
version of the shifted equations is 

g(t)^-;3/S- <" + ''"^ + ^' 7- , °f ^ , S + .Mt), (19a) 

a (a + /i) (7 + /i) 

m = f3IS + i^^±i^Kl±i^/ + f^[-f^ -i- + f^)il + f^)] s _ + + a^Mt), (19b) 

a (a + /i)(7 + /i) 

J(t) = «E-(7 + /i)J + (T303W. (19c) 

As Eqs. fll9ap - fll9cp are transformed using Eqs. (17al) - (17c|) . the a = ao//i scaling, the 
7 = 7o//U scaling, and the t = fir time scaling, the noise also is scaled so that the stochastic, 
transformed equations are given by 

dU /i^ ^^^Y _ ^^fj) + ^2) ^ ^2) [(^^ + 7o) f/ + -foW] 

— aQiJ + 



dr ao + 7o «o («o + 7o) 

^loW + («o + lo)U + , , {U + V) + a404, (20a) 



«o + 7o V (70 + A*^) («o + At^ 

dV ii"^ {-foV - aoU) (70 + /i^) («o + Ai^) [(ao + 7o) f/ + 7oW^] 

—aQU 



dr ao + 7o 7o ("o + 7o) 

'^loW + (ao + 70) f/ + , , , ) (f/ + V) + ^505, (20b) 



7o (ao + 7o) V (70 + At^) (ao + /i^) 

dr ao7o 



2ta , /^/^ TT/ , , - N rr , /^^ao7o 



/iV + ^ 70 + (ao + lo)U+ . ^ (f/ + ^) + ^606, (20c) 

70 V (70 + /^^) (ao + /i^) / 

where 

0-404 = q'202, (21a) 

ao + 7o 

£7505 =Aiai0i H ^^^^^0-202, (21b) 

ao + 7o 

0-606 =AiC^202 + AiO-303- (21c) 

The stochastic terms 04, 05, and 06 in Eqs. fl20ap - fl20cp are still additive, Gaussian noise 
processes. However, Eqs. ( I21ap - (l21cl) show how the transformation has acted on the original 
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stochastic terms 0i, 02, and 03 to create new noise processes which have a variance different 
from that of the original noise processes. Also note that we have suppressed the argument 
of 04, 05, and 06 in Eqs. f l20ap - fl20cl) . The time scaling means that these noise terms should 
be evaluated at /xr. 

X 10"* 




20 40 60 80 100 

t 



FIG. 2: (Color online) Time series of the fraction of the population that is infected with a disease, 
/, computed using the original, stochastic system of equations of the SEIR model [Eqs. (jlap -(fTc ] )] 
(red, solid line), and computed using the transformed, stochastic system of equations of the SEIR 
model [Eqs. (|20ap - (j20cp ] (blue, dashed line). The standard deviation of the noise intensity used in 
the simulation is (Tj = 0.0005, i = 1, . . . , 6. 

The system of equations given by Eqs. fl20ap - fl21cp are an exact transformation of the 
system of equations given by Eqs. f|Taj) - f|Tc|) . We numerically integrate the original, stochastic 
system of the SEIR model [Eqs. (fTa]) - (fTc|) ] along with the transformed, stochastic system 
[Eqs. f l20al) - fl20cp ] using a stochastic fourth-order Runge-Kutta scheme with a constant time 
step size. The original system is solved for S, E, and /, while the transformed system is 
solved for U, V, and W. In the latter case, once the values of U, V, and W are known, we 
compute the values of S, E, and / using the transformation given by Eqs. fl7aj) - fl7c|) . We 
shift S, E, and / respectively by 5*0, Eq, and Iq to find the values of S, E, and /. 

Figure [2] compares the time series of the fraction of the population that is infected with 
a disease, /, computed using the original, stochastic system of equations of the SEIR model 
with the time series of I computed using the transformed, stochastic system of equations of 
the SEIR model. 
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Although the two time series shown in Fig. [2] generally agree very well, there is some 
discrepancy. This discrepancy is due to the fact that the noise processes (7404, (J^cj)^, and 
(T606 of the transformed system are new, independent noise processes with different variance 
than the o"i</)i, cr202; and ^^(j)^ noise processes found in the original system. If we carefully 
take the noise realization from the original system, transform this noise using Eqs. fl21ap - 
fl21cp . and use this realization to solve the transformed system, then the two solutions would 
be identical. 



B. Reduction of the stochastic SEIR model 



It is tempting to consider the reduced stochastic model found by substitution of Eq. (fT7|l 

into Eqs. fl20bp and fl20cp . so that one has the following stochastic evolution equations (that 

hopefully describe the dynamics on the stochastic center manifold): 
_ _ /i^7o^aoVF _ fi^apW _ -fpfj.'^V _ (70 + /x^) Wap 
dr (ao + 7o)^ (oq + 7o)^ "0 + 7o "o + 7o 

(ao + 7o)' V (70 + /i^) («o + ^^')J\ {ao + 70) V 

dW fi^W , fi'W 

—r- =-, 72 ^ 

dr (ao + 7o) ao + 70 

( w + , + { V - + (22b) 



"0 + 70 V (70 + At^) ("0 + /W^) / V (ao + 7o)' 
One should note that Eqs. fl22al) and fl22bl) also can be found by naively adding the 

stochastic terms to the reduced system of evolution equations for the deterministic problem 
[Eqs. (11 Sal) and fllSbl) ]. This type of stochastic center manifold reduction has been done 
for the case of discrete noise^I. Additionally, in many other fields (e.g. oceanography, solid 
mechanics, fluid mechanics), researchers have performed stochastic model reduction using 
a Karhunen-Loeve expansion (principal component analysis, proper orthogonal decomposi- 
tion)^!^. However, this linear projection does not properly capture the nonlinear effects. 
Furthermore, one must subjectively choose the number of modes needed for the expansion. 
Therefore, even though the solution to the reduced model found using this technique may 
have the correct statistics, individual solution realizations will not agree with the original, 
complete solution. 

We will show that Eqs. fl22al) and (122bp do not contain the correct projection of the noise 
onto the center manifold. Therefore, when solving the reduced system, one does not obtain 
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the correct solution. Such errors in stochastic epidemic modehng impact the prediction of 
disease outbreak when modehng the spread of a disease in a population. 




FIG. 3: (Color online) Time series of the fraction of the population that is infected with a disease, 
/, computed using the complete, stochastic system of transformed equations of the SEIR model 
[Eqs. ()20ap - ()20cp ] (red, solid line), and computed using the reduced system of equations of the 
SEIR model that is based on the deterministic center manifold with a replacement of the noise 
terms [Eqs. (j22ap and ([22b])] (blue, dashed line). The standard deviation of the noise intensity 
used in the simulation is ai = 0.0005, i = 4,5, 6. The time series is shown for (a) i = to t = 40, 
and for (b) t = AO to t = 100. 

Using the same numerical scheme previously described, we numerically integrate the 
complete, stochastic system of transformed equations of the SEIR model [Eqs. fl20ap - fl20cl) 
along with the reduced system of equations that is based on the deterministic center manifold 
with a replacement of the noise terms [Eqs. ( \22a\] and (I22bp ]. The complete system is solved 
for U, V, and W, while the reduced system is solved for V and W. In the latter case, U is 
computed using the center manifold equation given by Eq. f|T7|) . Once the values of U, V, 
and W are known, we compute the values of S, E, and I using the transformation given by 



15 



Eqs. (17al) - (17c|) . We shift S, E, and I respectively by 5*0, Eq, and /q to find the values of S, 
E, and /. 

Figures [3]^a)-(b) compares the time series of the fraction of the population that is infected 
with a disease, /, computed using the complete, stochastic system of transformed equations 
of the SEIR model [Eqs. fl20al) - fl20cp ] with the time series of I computed using the reduced 
system of equations of the SEIR model that is based on the deterministic center manifold 
with a replacement of the noise terms [Eqs. fl22ap and fl22b|) ]. Figure [3]|^a) shows the initial 



transients, while Fig. [3]^b) shows the time series after the transients have decayed. One can 
see that the solution computed using the reduced system quickly becomes out of phase with 
the solution of the complete system. Use of this reduced system would lead to an incorrect 
prediction of the initial disease outbreak. Additionally, the predicted amplitude of the initial 
outbreak would be incorrect. The poor agreement, both in phase and amplitude, between the 
two solutions continues for long periods of time as seen in Fig. [3](b). We also have computed 
the cross-correlation of the two time series shown in Fig. [3]^a)-(b) to be approximately 0.34. 
Since the cross-correlation measures the similarity between the two time series, this low 
value quantitatively suggests poor agreement between the two solutions. 

Using the same systems of transformed equations, we compute 140 years worth of time 
series for 500 realizations. Ignoring the first 40 years of transient solution, the data is used 
to create a histogram representing the probability density, psi of the 5* and / values. Fig- 
ure |4](a) shows the histogram associated with the complete, stochastic system of transformed 
equations, while Fig. 111(b) shows the histogram associated with the reduced system of equa- 
tions with a replacement of the noise terms. The color-bar values in Figs. 111(a)- (b) have been 
normalized by 10 

One can see by comparing Fig. 111(a) with Fig. 111(b) that the two probability distributions 
qualitatively look the same. It also is possible to compare the two distributions using a 
quantitative measure. The KuUback-Leibler divergence, or relative entropy, measures the 
difference between the two probability distributions as 



(23) 



where Pij refers to the (i,j)th component of the probability density found using the com- 
plete, stochastic system of transformed equations [Fig. 111(a)], and Qij refers to the (i,j)th 
component of the probability density found using the reduced system of equations [Fig. Hfb)] . 

16 



X 10 




0.058 0.062 0.066 0.07 
S 

X 10""^ 




0.058 0.062 0.066 0.07 

S 

FIG. 4: (Color online) Histogram of probability density, psi of the S and / values found using (a) 
the complete, stochastic system of transformed equations for the SEIR model [Eqs. (j20ap - (j20cp ]. 
and (b) the reduced system of equations of the SEIR model that is based on the deterministic 
center manifold with a replacement of the noise terms [Eqs. (j22ap and (|22bp ]. The histograms are 
created using 100 years worth of time series (starting with year 40) for 500 realizations, and the 
color-bar values have been normalized by 10~^. 

In our numerical computation of the relative entropy, we have added 10~^° to each Pij and 
Qij. This eliminates the possibility of having a Qij = in the denominator of Eq. (123!) and 
does not have much of an effect on the relative entropy. 

If the two histograms were identical, then the relative entropy given by Eq. (125]) would 
he dxL = 0. The two histograms shown in Figs. |4](a)-(b) have a relative entropy oi dxL = 
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0.0391, which means that the two histograms, while not identical, are quantitatively very 
similar. However, one cannot rely entirely on the histograms alone to say that the solutions 
of the complete system and the reduced system agree. As we have seen in Figs. IS]^a)-(b), 
the two solutions have differing amplitudes and are out of phase with one another. It is 
important to note that these features are not picked up by the histograms of Fig. |H 



V. CORRECT PROJECTION OF THE NOISE ONTO THE STOCHASTIC CEN- 
TER MANIFOLD 

To project the noise correctly onto the center manifold, we will derive a normal form 
coordinate transform for the complete, stochastic system of transformed equations of the 
SEIR model given by Eqs. fl20ap - fl20cl) . The particular method we use to construct the 
normal form coordinate transform not only reduces the dimension of the dynamics, but also 
separates all of the fast processes from all of the slow processes^. This technique has been 
modified and applied to the large fluctuations of multiscale problems^^. 

Many publicationsi^'^>2ii2^ discuss the simplification of a stochastic dynamical system 
using a stochastic normal form transformation. In some of this work— i^^, anticipative noise 
processes appeared in the normal form transformations, but these integrals of the noise 
process into the future were not dealt with rigorously. 

Later, the rigorous, theoretical analysis needed to support normal form coordinate trans- 



forms was developed in Refs. 
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J24l . The technical problem of the anticipative noise integrals 
also was dealt with rigorously in this work. Even later, another stochastic normal form trans- 
formation was developed^^. This new method allows for the "[removal of] anticipation ... 
from the slow modes with the result that no anticipation is required after the fast transients 



decay" (Ref. l25|, pp. 13). The removal of anticipation leads to a simplification of the normal 



form. Nonetheless, this simpler normal form retains its accuracy with the original stochastic 
system^^. 

We shall use the method of Ref. 
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to simplify our stochastic dynamical system to one 
that emulates the long-term dynamics of the original system. The method involves five 
principles, which we recapitulate here for completeness: 

1. Avoid unbounded, secular terms in both the transformation and the evolution equa- 
tions to ensure a uniform asymptotic approximation. 
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2. Decouple all of the slow processes from the fast processes to ensure a valid long-term 
model. 



3. Insist that the stochastic slow manifold is precisely the transformed fast processes 
coordinate being equal to zero. 

4. To simplify matters, eliminate as many as possible of the terms in the evolution equa- 
tions. 

5. Try to remove all fast processes from the slow processes by avoiding as much as possible 
the fast time memory integrals in the evolution equations. 

In practice, the original stochastic system of equations (which satisfies the necessary spec- 
tral requirements) in (f/, V, W)'^ coordinates is transformed to a new {Y, Xi, X2)^ coordinate 
system using a near-identity stochastic coordinate transform given as 



where the specific form of ^ (F, Xi, X2, r), rj (F, Xi, X2, r), and p (F, Xi, X2, r) is chosen to 
simplify the original system according to the five principles listed previously, and is found 
using an iterative procedure. To outline the procedure, we provide details for a simple 
example in Appendix 1X1 

Several iterations lead to coordinate transforms for f/, V , and W along with evolution 
equations describing the F-dynamics, Xi-dynamics, and X2-dynamics in the new coordinate 
system. The F-dynamics have exponential decay to the F = slow manifold. Substitution 




(24a) 




(24b) 



Py=X2 + p(F,Xi,X2,r), 



(24c) 
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of y = leads to the coordinate transforms 
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All of the stochastic terms in Eqs. f l25ap -( l25cl) consist of integrals of the noise process 
into the past (convolutions), as given by Eqs. (126!) and (!27|) . These memory integrals are 
fast-time processes. Since we are interested in the long-term slow processes and since the 
expectation of Q equals e"^"^ * E[(f)\, where -^[0] = 0, we neglect the memory integrals and 
the higher-order multiplicative terms found in Eqs. fl25ap - fl25cl) so that 
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FIG. 5: (Color online) Time series of the fraction of the population that is infected with a disease, 
/, computed using the complete, stochastic system of transformed equations of the SEIR model 
[Eqs. (j20ap - (j20c|) ] (red, solid line), and computed using the reduced system of equations of the 
SEIR model that is found using the stochastic normal form coordinate transform [Eqs. (j29ap . 
()29bp . (jBlap . and (|Blbp ] (blue, dashed line). The standard deviation of the noise intensity used in 
the simulation is dj = 0.0005, i = 4,5, 6. The time series is shown for (a) t = to t = 40, and for 
(b) t = 40 to t = 100. 



Note that Eq. fl28ap is the deterministic center manifold equation, and at first-order, matches 
the center manifold equation that was found previously [Eq. (IT7|) ]. 

Substitution of F = and neglecting all multiplicative noise terms and memory integrals 
using the argument from above (so that we consider only first-order noise terms) leads to 
the following reduced system of evolution equations on the center manifold: 

dXi 



dr 
dX2 

dr 



F(Xi(r),X2(r)), 
G(Xi(r),X2(r)). 



(29a) 
(29b) 



The specific form of F and G in Eqs. fl29al) and (129bl) are complicated, and are therefore 
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presented in Appendix [Bl 

We numerically integrate the complete, stochastic system of transformed equations of the 
SEIR model [Eqs. fl20al) - fl20cp ] along with the reduced system of equations that is found using 
the stochastic normal form coordinate transform [Eqs. fl29ap . fl29bp . fIBlap . and flBlb|) ]. The 
complete system is solved for U, V, and W, while the reduced system is solved for Xi = V 
and X2 = W . In the latter case, U is computed using the center manifold equation given 
by Eq. fl28al) . As before, once the values of f/, V , and W are known, we compute the values 
of S, E , and / using the transformation given by Eqs. (17al) - (!7cl) . We shift S, E, and I 
respectively by Sq, Eq, and Jq to find the values of S, E, and /. 

Figures [ni^a)-(b) compares the time series of the fraction of the population that is infected 
with a disease, /, computed using the complete, stochastic system of transformed equations 
of the SEIR model [Eqs. fl20al) - fl20cp ] with the time series of / computed using the reduced 
system of equations of the SEIR model that is found using the stochastic normal form 
coordinate transform [Eqs. ( 129ap . ( ]29b[) . ( IBlaj) . and ( IBlbl) ]. Figure [5]^a) shows the initial 
transients, while Fig. Mjo) shows the time series after the transients have decayed. One 
can see that there is excellent agreement between the two solutions. The initial outbreak is 
successfully captured by the reduced system. Furthermore, Fig. [5](b) shows that the reduced 
system accurately predicts recurrent outbreaks for a time scale that is orders of magnitude 
longer than the relaxation time. This is not surprising since the solution decays exponentially 
throughout the transient and then remains close to the lower- dimensional center manifold. 
Since we are not looking at periodic orbits, there are no secular terms in the asymptotic 
expansion, and the result is valid for all time. Additionally, any noise drift on the center 
manifold results in bounded solutions due to sufficient dissipation transverse to the manifold. 
The cross-correlation of the two time series shown in Fig. [5] is approximately 0.98, which 
quantitatively suggests there is excellent agreement between the two solutions. 

Using the same systems of transformed equations, we compute 140 years worth of time 
series for 500 realizations. As before, we ignore the first 40 years worth of transient solution, 
and the data is used to create a histogram representing the probability density, psi of the 
S and / values. Figure [6](a) shows the histogram associated with the complete, stochastic 
system of transformed equations, while Fig. [6](b) shows the histogram associated with the 
reduced system of equations found using the normal form coordinate transform. The color- 
bar values in Figs. [6](a)-(b) have been normalized by 10~^. 
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FIG. 6: (Color online) Histogram of probability density, "psi of the S and / values found using 
(a) the complete, stochastic system of transformed equations for the SEIR model with mortality 
[Eqs. (|2Uap - (j20cp ]. and (b) the reduced system of equations of the SEIR model with mortality 
that is found using the stochastic normal form coordinate transform [Eqs. (j29ap . (|29bp . (|Blap . 
and (jBlbp ]. The histograms are created using 100 years worth of time series (starting with year 
40) for 500 realizations, and the color-bar values have been normalized by 10"'^. 

As we saw with Figs. 111(a)- (b), the probability distribution shown in Fig. [n](a) looks 
qualitatively the same as the probability distribution shown in Fig.[6](b). Using the Kullback- 
Leibler divergence given by Eq. (l23ll . we have found that the two histograms shown in 
Figs. [6]|^a)-(b) have a relative entropy of (Ikl = 0.0953. Since this value is close to zero, the 
two histograms are quantitatively very similar. 
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FIG. 7: Cross-correlation between time series found using the original, stochastic system of trans- 
formed equations and the reduced system of equations based on the deterministic center manifold 
("circle" markers), and cross-correlation between time series found using the original, stochastic 
system of transformed equations and the reduced system of equations based on the stochastic nor- 
mal form coordinate transform ("square" markers). The cross-correlation is computed using time 
series from t = 800 to t = 1000. 

In addition to computing the cross-correlation between the solution of the original system 
and the solutions of the two reduced systems for cTj = 0.0005, we have computed the cross- 
correlation for the case of zero noise as well as for noise intensities ranging from a = 5.0 x 
10~^° to cr = 5.0 X 10~^. These cross-correlations were computed using time series from 
t = 800 to t = 1000. For the deterministic case (zero noise), the cross-correlation between 
the time series which were computed using the original system and the reduced system 
based on the deterministic center manifold is 1.0, since the agreement is perfect. The cross- 
correlation between the original system and the reduced system found using the stochastic 
normal form is also 1.0. Figure [7] shows the cross-correlation between the original system 
and the two reduced systems for various values of a. 

One can see in Fig. [7] that the solutions found using the reduced system based on the 
deterministic center manifold compare poorly with the original system at very low noise 
values. Furthermore, as the noise increases, the agreement between the two solutions gets 
worse. On the other hand, Fig.[7]shows that the solutions computed using the reduced system 
found using the normal form coordinate transform compare very well with the solutions to 
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the original system across a wide range of small noise intensities. 
VI. DISCUSSION 

We have demonstrated that the normal form coordinate transform method reduces the 
Langevin system so that both the noise and dynamics are accurately projected onto the 
lower-dimensional manifold. It is natural to consider (a) the replacement of the stochastic 
term by a deterministic, period drive of small amplitude, and (b) the extension to finite 
populations. These cases are discussed respectively in Sec. IVI Al and Sec. IVIBI 

A. The Case of Deterministic Forcing 

A single time series realization of the noise might be thought of as a deterministic function 
of small amplitude driving the system. One could rederive the normal form for such a 
deterministic function. However, since our derived normal form holds specifically for the 
case of white noise, we show that a simple replacement of the stochastic realization with a 
deterministic realization does not work. As an example, one could consider the following 
sinusoidal functions: 

a^cj)^ = cos (lOvr/it) /8000, (30a) 
a202 = sin {An fit) /8000, (30b) 
(T3(^3 = cos (lOvr/it) /8000, (30c) 

where (7404, cr^cp^, and ae^e are given by Eqs. fl21ap - fl21cp . Using Eqs. fl30al) - fl30cl) or some 
other similar deterministic drive, the solution computed using the reduced system based on 
the deterministic center manifold analysis will agree perfectly with the solution computed 
using the complete system of equations. On the other hand, since the reduced system based 
on the normal form analysis was derived specifically for white noise, the transient solution 
found using this reduced system will not agree with the solution found using the complete 
system. It is possible to find a normal form coordinate transform for periodic forcing, but 
the normal form will be different than the one derived in this article for white noise. 

Figures [8t^a)-(b) compares the time series of the fraction of the population that is infected 
with a disease, J, computed using the complete system of transformed equations of the 
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FIG. 8: (Color online) Time series of the fraction of the population that is infected with a disease, 
/, computed using the complete system of transformed equations of the SEIR model [Eqs. (j20ap - 
(j20cp ] (red, solid line), and computed using the reduced system of equations of the SEIR model 
that is found using the normal form coor dinate transform [Eqs. <^M), (j29b]) . ([BTa]) . and ([Bib]) ] 
(blue, dashed line). The stochastic terms in both systems have been replaced by the deterministic 
terms given by Eqs. (|30ap - ()30bp . The time series is shown from (a) t = to t = 25, and from (b) 
t = 65 to t = 70. 



SEIR model [Eqs. (120ap - fl20cp ] with the time series of / computed using the reduced system 
of equations of the SEIR model that is found using the stochastic normal form coordinate 
transform [Eqs. fl29ap . (]29bp . (IB lap , and (IBlbp ], but where the stochastic terms of both 
systems have been replaced by the deterministic terms given by Eqs. (I30al) - (l30cl) . Figure [Hl^a) 
shows the initial transients, while Fig. [8]^b) shows a piece of the time series after the transients 
have decayed. One can see in Figs. [8]^a)-(b) that although the two solutions eventually 
become relatively synchronized with one another, there is poor agreement, both in phase 
and amplitude, throughout the transient. 
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B. The Case of Finite Populations 



The solutions to the original system and both reduced systems are continuous solutions 
based on an infinite population assumption, and are found using Langevin equations having 
Gaussian noise. It is interesting to examine the effects of general noise by using a Markov 
simulation to compare solutions of the original and reduced systems. 

The complete system in the original variables (see page 5) will evolve in time t in the 
following way: 



transition 


rate 
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(31) 



Using a total population size of = 10 million, we have performed a Markov simulation of 
the system. After completing the Markov simulation, we divided s, e, and i by to find S, 
E, and I. Figure [9]^a) shows a time series, after the transients have decayed, of the fraction 
of the population that is infected with a disease, I. The results refiect both the mean and 
the frequency of the deterministic system. Performing the simulation for 500 realizations 
allows us to create a histogram representing the probability density, psi of the 5* and / 
values. This histogram is shown in Fig. Mjo), and one can see that the probability density 
refiects the amplitude, which varies with the population size, of S and /. The color-bar 
values in Fig. [9](b) have been normalized by 10~^ 

The complete system in the transformed variables has the stable endemic equilibrium at 
the origin. To bound the dynamics to the first octant, we use the fact that s > 0, e > 0, and 
i > to derive the appropriate inequalities for the transformed, discrete variables u, v, and 
w. These inequalities can be found in Appendix [C] as Eq. flCl|) . These inequalities enable us 
to define new discrete variables Yi, Y2, and Y3, given by Eqs. (lC2ap - (lC2cP in Appendix O 

In the Yi variables, we define evolution relationships similar to those found in Eq. flMl) . 
The complete transformed system will evolve in time r according to the transition and rates 
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FIG. 9: (a) Time series of the fraction of the population that is infected with a disease, /, computed 
using a Markov simulation of the complete, original equations of the SEIR model [Eq. pip ], and (b) 
(color online) a histogram of probability density, psi of the S and / values found using a Markov 
simulation of Eq. (j3ip . The histogram is created using 100 years worth of data (starting with year 
40) for 500 realizations, and the color-bar values have been normalized by 10~^. 

given by Eq. (1C3P in Appendix O 

After performing a Markov simulation of Eq. (IC3p with a population size of = 10 
million, we can compare the dynamics of the transformed system to the dynamics of the 
original system by transforming the Yi variables in the time series back to the original s, 
e, and i variables. Dividing by A^ yields S*, and I. Figure [TO](a) shows a time series, 
after the transients have decayed, of the fraction of the population that is infected with a 
disease, /. The mean and the frequency agree with those found from the Markov simula- 
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FIG. 10: (a) Time series of the fraction of the population that is infected with a disease, I, 
computed using a Markov simulation of the complete, transformed equations of the SEIR model 
[Eq. (|C3p ]. and (b) (color online) a histogram of probability density, psi of the S and / values 
found using a Markov simulation of Eq. ()C3p . The histogram is created using 100 years worth of 
data (starting with year 40) for 500 realizations, and the color-bar values have been normalized by 
10-1 

tion of the original system. We have performed the simulation for 500 realizations, and a 
histogram representing the probability density, psi is shown in Fig. [TU](b). The color-bar 
values in Fig. [TUT b) have been normalized by 10"^. One can see in Fig. [TU](a) that the 
relative fluctuations of the / component has nearly doubled. While the fluctuation size was 
0.152 for the original system, it is 0.310 for the transformed system. Additionally, the two 
histograms shown in Figs. [9](b) and [TOVb) have a relative entropy of (Ikl = 0.9519, which 
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means they are not in agreement. Because the simulation of the stochastic dynamics in the 
complete system of transformed variables do not qualitatively (or quantitatively) resemble 
the original stochastic system, we cannot expect that the reduced system will agree with 
either the original or the transformed systems. Therefore, much care should be exercised 
when extending the model reduction results (which show outstanding agreement) derived 
for a specific type of noise in the limit of infinite population to finite populations with a 
more general type of noise. 

VII. CONCLUSIONS 

We have considered the dynamics of an SEIR epidemiological model with stochastic 
forcing in the form of additive, Gaussian noise. We have presented two methods of model 
reduction, whereby the goal is to project both the noise and the dynamics onto the stochastic 
center manifold. The first method uses the deterministic center manifold found by neglecting 
the stochastic terms in the governing equations, while the second method uses a stochastic 
normal form coordinate transform. 

Since the original system of governing equations does not have the necessary spectral 
structure to employ either deterministic or stochastic center manifold theory, the system 
of equations has been transformed using an appropriate linear transformation coupled with 
appropriate parameter scaling. At this stage, the first method of model reduction can be 
performed by computing the deterministic center manifold equation. Substitution of this 
equation into the complete, stochastic system of transformed equations leads to a reduced 
system of stochastic evolution equations. 

The solutions of the complete, stochastic system of transformed equations as well as 
the reduced system of equations were computed numerically. We have shown that the 
individual time series do not agree, because the noise has not been correctly projected onto 
the stochastic center manifold. However, by comparing histograms of the probability density, 
Psi of the S and / values, we saw that there was very good agreement. This is caused by 
the fact that although the two solutions are out of phase with one another, their range of 
amplitude values are similar. The phase difference is not represented in the two histograms. 
This is a real drawback when trying to predict the timing of outbreaks, and leads to potential 
problems when considering epidemic control, such as the enhancement of disease extinction 
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through random vaccine control^. 

To accurately project the noise onto the manifold, we derived a stochastic normal form 
coordinate transform for the complete, stochastic system of transformed equations. The 
numerical solution to this reduced system was compared with the solution to the original 
system, and we showed that there was excellent agreement, both qualitatively and quanti- 
tatively. As with the first method, the histograms of the probability density, psi oi the 5* 
and / values agree very well. 

It should be noted that the use of these two reduction methods is not constrained to 
problems in epidemiology, but rather may be used for many types of physical problems. For 
some generic systems, such as the singularly perturbed, damped Duffing oscillator, either 
reduction method can be used since the terms in the normal form coordinate transform which 
lead to the average stochastic center manifold being different from the deterministic center 
manifold occur at very high order—. In other words, the average stochastic center manifold 
and deterministic center manifold are virtually identical. For the SEIR model considered 
in this article, there are terms at low order in the normal form transform which cause a 
significant difference between the average stochastic center manifold and the deterministic 
manifold. Therefore, as we have demonstrated, when working with the SEIR model, one 
must use the normal form coordinate transform method to correctly project the noise onto 
the center manifold. 

In summary, we have presented a new method of stochastic model reduction that allows 
for impressive improvement in time series prediction. The reduced model captures both 
the amplitude and phase accurately for a temporal scale that is many orders of magnitude 
longer than the typical relaxation time. Since sufficient statistics of disease data are limited 
due to short time series collection, the results presented here provide a potential method to 
properly model real, stochastic disease data in the time domain. Such long-term accuracy of 
the reduced model will allow for the application of effective control of a disease where phase 
differences between outbreak times and vaccine controls are important. Additionally, since 
our method is general, it may be applied to very high-dimensional epidemic models, such 
as those involving adaptive networks. From a dynamical systems viewpoint, the reduction 
method has the potential to accurately capture new, emergent dynamics as we increase 
the size of the random fluctuations. This could be a means to identify new noise-induced 
phenomena in generic stochastic systems. 
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APPENDIX A: DETAILS OF THE ITERATIVE PROCEDURE FOR A SIMPLE 
EXAMPLE 

We consider the system given by 

dx 

— =/i(y + a0), (Ala) 

^=x-x^-y, (Alb) 
dr 

^ =0. (Ale) 
dr 

The iterative procedure begins by letting 

x^X, (A2a) 

X' ^ 0, (A2b) 
and by finding a change to the y coordinate (fast process) with the form 

y = Y + r]{T,X,Y) + ..., (A3a) 

Y' = -Y + G{t,X,Y) + (A3b) 

where rj and G are small corrections to the coordinate transform and the corresponding 
evolution equation. Substitution of Eqs. (]A2aP - (]A3bl) into Eq. f lAlbO gives the equation 



dr] dr] dX dr] dY , . ,^ 
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Replacing Y' = dY/dr with -Y + G [Eq. flMbl) ]. noting that dX/dr = [Eq. flMbl) ]. and 
ignoring the term drj/dY ■ G since it is a product of small corrections leads to 

Equation (lASp must now be solved for G and 77. In order to keep the evolution equation 
[Eq. (1A3bp ] as simple as possible (principle (j4]) of Sec. |V]), we let G = 0, which means 
that the coordinate transform [Eq. flASal) ] is modified hj rj = X — X^. Therefore, the new 
approximation of the coordinate transform and its dynamics are given by 

y = Y + X-X^ + 0{C^), (A6a) 

Y' = -Y + 0{C^), (A6b) 

where ( = \ {X, Y, fi, a)\ so that ( provides a count of the number of X, Y, fi, and a factors 
in any one term. 

For the second iteration, we seek a correction to the x coordinate (slow process) with the 
form 

x = X + ^{t,X,Y) + ..., (A7a) 
X' = F{t,X,Y) + (A7b) 

where ^ and F are small corrections. Substitution of Eqs. (IA6aP - (IA7bp into Eq. (lAlaP leads 

to 

Replacing X' = dX/dr with F [Eq. flATbl) ]. replacing dY/dr with -Y [Eq. flA6bl) ]. and 
ignoring the term d^/dX ■ F since it is a product of small corrections gives the equation: 

F+^-Y^ = ^x{Y + X-X^)+ ixa(t>. (A9) 

Equation (1A9P must now be solved for F and ^. As in the first step, we employ princi- 
ple dl]) and keep the evolution equation [Eq. (lA7bp ] as simple as possible. However, since 
the terms ^{X — X^) located on the right-hand side of Eq. flA9P do not contain r or F, these 
terms must be included in F. Therefore, one piece of F will be F = yu(X — X^). 

The remaining deterministic term on the right-hand side of Eq. ( 1A9P contains Y . This 
term can therefore be integrated into ^. The equation to be solved is 

-Y-^ = fiY, (AlO) 
33 



whose solution is given as ^ = ~jjY. 

To abide by principle (jl]), we would like to integrate the stochastic piece on the right-hand 
side of Eq. (lAOp into ^, by solving the equation 

d^/dr = fia(f). (All) 

However, the solution of Eq. (lAlip is given by 

^ = /i(T / 0rfr, (A12) 



which has secular growth like a Wiener process. Since this would violate principle ([T]), we 
must let F = ftacp. 

Putting the three pieces together yields ^ = —fiY and F = fi{X — X^) + fia(j). Therefore, 
the new approximation of the coordinate transform and its dynamics are given by 

X = X - fiY + 0{C^), (Al3a) 

X' = ij{X -X^)+fia(j)+0{C). (Al3b) 

The construction of the normal form continues by seeking corrections, ^ and F, to the 
X coordinate transform and the X evolution using the updated residual of the x equation 
[Eq. ( lAlaP ], and by seeking corrections, rj and G, to the y coordinate transform and the Y 
evolution equation using the updated residual of the y equation [Eq. (lAlbj) . 

APPENDIX B: REDUCED, STOCHASTIC SEIR MODEL: CORRECT PROJEC- 
TION OF THE NOISE 

The specific form of F and G in Eqs. f l29al) and (]29bp are given as 

F = - \ah'oX2 + r^^^l'o ( + c^oX.X,) + /i^ (ao7o'^i + 

(ao + 7o) V "o + 7o / 

«o7o (2ao + 5ao7o + 5ao7o + 7o) ^ «o/^'7o ^2 y «o/^'7o ^ , 

A2 — -, -A.^yi2 ~ ~ Tl' 1 2 I + 



(ao + 7of («o + 7o) (ao + 7o) 

aoPlo {a'o - ao7o - 3ao7o " 37;^) 



("0 + 7of 



I [70 («0 + 7o) ("0 + yU^) (70 + yU^)] + 



/i^(ao + "o7o + 7o) Z' ^40004 , cr67o0( 

Cr505 



(tto + 7o)^ V 7o "0 + 7o 

At^«o/3 / q"4ao04 ^ O'67O06 

("0 + 7o)^ V 7o ("0 + 7o) 



(Bla) 
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G 



.("o + 7o) (ao + 7o)^ / V ("o + 7o) 



(«o + 7o) (ao + 7o)^ / V (ao + 7o) 

7 ^ ^ Xi - -, ■ 7X2 - 3aoP7o («o + ) (7o + ) X1X2+ 

(tto + 7o) ("0 + 7o) 

al(3lo {al + 2ao7o + 



("0 + 7o) 



-XiX'i 



I [ao7o ("0 + yu^) (70 + y"^)] + o-6</'6+ 



(a^ + ao7o + 7o) , , erg [7^ (ao + /i^) (70 + At ) + ao7o («o + 7o)] , , 

1 \ N H ; 73 </>6+ 

ao7o (ao + 7o) ao (ao + 7o) 

/ (T404 , O-67O06 



(Bib) 



(q;o + 7o) V 70 (ao + 7o)' 
APPENDIX C: MARKOV SIMULATION FOR TRANSFORMED SEIR MODEL 

The complete system in the transformed variables has the stable endemic equilibrium at 

the origin. To bound the dynamics to the first octant, we transform the new variables by 

using the original properties of s > 0, e > 0, and z > 0, so that 

/i^A^7o A^7o + ao^ + 7o«o) , X (ao + 70) , x 

— -, — - — r, IT-, — - — ^ (CI) 

«o ("0 + 7o ) "oPAi («o + 7o ) 70^0 

Therefore, we define the following new variables: 

y, = -„ + ^^. (C2a) 

«o («o + 7o) 

V - „ ^ ^70 (/^/^^ + Qq^ + 7oao) /p^, x 

aoPn (ao + 7o) 

70^0 

In these variables, we define evolution relationships similar to Eq. (l3Ti) . The complete 
transformed system will evolve in r the following way: 

transition rate 

(^1 + 1.^2,13) ^{^y^Ys + Y,') 

(Fi - 1, Y2, F3) K + /i')n + f (^>^i>^3 + Y.Y^) 

(1^,^2 + 1,^3) /i^AT + f (^FiF3 + ^VIF^) . (C3) 

(Fi, - 1, 1^3) «on + f^^Y, + f ((^^2^3 + ^Y,' 

(n, F2, 13 + 1) («o + 70)1^1 + f (1^213 + ^^Fi^) 

(Fi, F2, ^3 - 1) (70 + fi')Ys + f (FiFs + ^^FiF^) 
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